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ABSTRACT 

We propose a new fast algorithm for solving one of the standard for- 
mulations of frame-based image deconvolution: an unconstrained 
optimization problem, involving an £2 data-fidelity term and a non- 
smooth regularizes Our approach is based on using variable splitting 
to obtain an equivalent constrained optimization formulation, which 
is then addressed with an augmented Lagrangian method. The result- 
ing algorithm efficiently uses a regularized version of the Hessian of 
the data fidelity term, thus exploits second order information. Exper- 
iments on a set of image deblurring benchmark problems show that 
our algorithm is clearly faster than previous state-of-the-art methods. 

1. INTRODUCTION 

1.1. Problem Formulation 

The standard model in image deblurring assumes that the noisy blurred 
observed version y, of an original image x £ R n , was obtained via 

y = Hx + w, 

where H is the matrix representation of a convolution and w is Gaus- 
sian white noise. In frame-based deblurring/deconvolution, the un- 
known image x is expressed as x = W/3, where the columns of 
matrix W are the elements of a frame, such as a wavelet orthonor- 
mal basis or a redundant dictionary [6], [7], [8], [11], [13], [14]. The 
coefficients of this representation are then estimated, under one of 
the well-known sparsity inducing regularizers, typically the l\ norm, 
leading to the optimization problem 

3 = argmi n i||HW/3-y||^+T0(/3); (D 
P z 

in (1), (j) '■ R m — > K is the regularizer, which is usually convex 
but nonsmooth, and r > is the regularization parameter [7]. This 
formulation is called the synthesis approach [12], since it is based 
on the synthesis equation x = W/3. In the last decade, a consid- 
erable amount of research has been devoted to designing efficient 
algorithms for solving (1). This interest has been further stimulated 
by the recent emergence of compressive sensing (CS) [5], [9], since 
CS reconstruction can be formulated as (1) [17], [24]. 

1.2. Previous Algorithms 

In most practical problems (including CS), matrix HW (and even 
H or W) cannot be stored explicitly and it is highly impractical 
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to access portions (lines/columns or blocks) of it. These facts pre- 
clude most off-the-shelf optimization algorithms from being directly 
used and has stimulated the development of special purpose meth- 
ods. These methods operate under the constraint H and W (and 
their transposes) can only be used to form matrix-vector products, 
since these products can be performed efficiently using the FFT and 
fast wavelet transforms. 

Arguably, the standard algorithm for solving (1) is the so-called 
iterative shrinkage/thresholding (1ST), which can be derived from 
different viewpoints: expectation-maximization [13], majorization- 
minimization [8], [14], forward-backward operator splitting [7], [16]. 
A key ingredient of 1ST is the so-called shrinkage/thresholding func- 
tion associated to </>, <J> T ^ : R n — > R n , defined as 

* T0 (a) = argmini||a - pf 2 + r<t>{0). (2) 
P z - 

An excellent coverage of these functions, also known as Moreau 
proximal maps, can be found in [7]. 

The fact that 1ST tends to be slow, in particular when H is poorly 
conditioned, has stimulated some recent research aimed at obtaining 
faster variants. The recent two-step 1ST (TwIST) algorithm [3], in 
which each iteration uses the two previous iterates (rather than only 
the previous one, as in 1ST), was shown to be considerably faster 
than 1ST on various deconvolution problems. Another two-step vari- 
ant of 1ST, named fast 1ST algorithm (FISTA), was recently proposed 
and also shown to be faster than 1ST [2]. A recent strategy to obtain- 
ing faster variants of 1ST consists in using more aggressive choices 
of step size in each iteration. This is the case in the SpaRSA (sparse 
reconstruction by separable approximation) framework [21], [22], 
which was also shown to clearly outperform standard 1ST. 

1.3. Proposed Approach 

The approach proposed in this paper is based on variable splitting. 
The idea is to split the variable /3 into a pair of variables (3 and 8, 
each to serve as the argument of each of the two functions in (1), and 
then minimize the sum of the two functions under the constraint that 
the two variables have to be equal, thus making the problems equiv- 
alent. This rationale has been recently used in the split-Bregman 
method [15], which was proposed to address constrained optimiza- 
tion formulations for solving inverse problems. In this paper, we 
exploit a different splitting to attack problem (1), arguably the most 
classical formulation for frame-based regularization of linear inverse 
problems [6], [7]. 

The constrained optimization problem produced by the splitting 
procedure is addressed using an augmented Lagrangian (AL) algo- 
rithm [18]. AL was shown to be equivalent to the Bregman itera- 
tive methods [19], [23]. We adopt the AL perspective, rather than 



the Bregman view, as it is a more standard optimization tool. We 
show that by exploiting the fact that W is a frame, the resulting 
algorithm solves (1) much faster than the previous state-of-the-art 
methods FISTA [2], TwIST [3], and SpaRSA [22]. 

The speed of the proposed algorithm may be justified by the fact 
that it uses (a regularized version of) the Hessian of the data fidelity 
term, W T H T H W, while the above mentioned algorithms essen- 
tially only use gradient information. Although, as referred earlier, 
this matrix can not be formed, we show that if W is a tight frame 
and H a convolution, our algorithm can use it in an efficient way. 

2. BASIC TOOLS 
2.1. Variable Splitting 

Consider an unconstrained optimization problem in which the objec- 
tive is the sum of two functions: 



min /i(u) + /2(u). 



(3) 



Variable splitting (VS) is a simple procedure in which a new variable 
v is introduced to serve as the argument of /2, under the constraint 
that u = v. In other words, the constrained problem 



mm 
s.t. 



/i(u) + / 2 (v) 



(4) 



v}, the 



is equivalent to (3), since in the feasible set {(u, v) : 
objective function in (4) coincides with that in (3). 

VS was used in [20] to derive a fast algorithm for total-variation 
based restoration. VS was also used in [4] to handle problems where 
instead of the single regularizer r0(/3) in (1), there is a linear combi- 
nation of two (or more) regularizes: ti0i(/3) + t~202(/3). In [4] and 
[20], the constrained problem (4) is attacked by a quadratic penalty 
approach, i.e., by solving 

min /i(u)+/ 2 (v) + g||u-v||i, (5) 

by alternating minimization with respect to u and v, while slowly 
increasing fi to force the solution of (5) to approach that of (4). The 
idea is that each step of this alternating minimization may be much 
easier than the original unconstrained problem (3). The drawback is 
that as /i increases, the intermediate minimization problems become 
increasingly ill-conditioned, thus causing numerical problems [18]. 

A similar VS approach underlies the recently proposed split- 
Bregman methods [15]. In those methods, the constrained problem 
(4) is addressed using a Bregman iterative algorithm, which has been 
shown to be equivalent to the AL method [23]. 

2.2. Augmented Lagrangian 

Consider a linear equality constrained optimization problem 



min E(z) 
s.t. Az - b =0, 



(6) 



where b £ W and A £ W xd . The so-called augmented La- 
grangian function for this problem is defined as 

£a(z, A, M ) = E(z) + A T (Az - b) + | || Az - b|||, (7) 

where A £ W is a vector of Lagrange multipliers and p, > is 
called the AL penalty parameter [18]. The AL algorithm iterates 



between minimizing £a(z, A, jx) with respect to z, keeping A fixed, 
and updating A. 

Algorithm AL 

1. Set k = 0, choose fi > 0, zo, and Ao. 

2. repeat 

3. z fc+ i 6argmin z £A(z,At,ji) 

4. A fc+ i <- A fc + (i(Az H1 - b) 

5. k^k + 1 

6. until stopping criterion is satisfied. 

It is possible (in some cases recommended) to update the value 
of /i at each iteration [18], [1] (Chap. 9). However, unlike in the 
quadratic penalty method, it is not necessary to take /i to infinity to 
guarantee that the AL converges to the solution of the constrained 
problem (6). In this paper, we will consider only the case of fixed p. 

After a straightforward manipulation, the terms added to E(z) 
in Ca(z, Afc, /i) (see (7)) can be written as a single quadratic term, 
leading to the following alternative form for the AL algorithm: 

Algorithm AL ( version 2 ) 

1. Set k — 0, choose /i > 0, zo, and do. 

2. repeat 

3. z fc+ i € argmin z £;(z) + f ||Az - d fc ||| 

4. d fe+ i <— d fc - (Az fc+ i — b) 

5. k^k + 1 

6. until stopping criterion is satisfied. 

This form of the AL algorithm makes clear its equivalence with 
the Bregman iterative method, as given in [23]. 

2.3. AL for Variable Splitting and Its Convergence 

Problem (4) can be written in the form (6) with z = [u T , v T ] T , 
b = 0, A = [I - I], and E(z) = /i(u) + / 2 (v). With these 
definitions in place, Steps 3 and 4 of the AL algorithm (version 2) 
can be written as follows 



-v -d fe ||| (8) 



( vT + \ ) ear s™ n /i( u ) + M v ) + fll u - v 

d fc+1 = d fc - (ufc+i - Vfc+i). (9) 
The minimization problem (8) is clearly non-trivial: in general, it 
involves non-separable quadratic and possibly non-smooth terms. A 
natural approach is to use a non-linear block-Gauss-Seidel (NLBGS) 
technique, in which (8) is solved by alternating minimization with 
respect to u and v, while keeping the other variable fixed. Remark- 
ably, it has been shown that the AL algorithm converges, even if 
the exact solution of (8) is replaced with a single NLBGS step [10, 
Theorem 8] (see also [19]). The resulting algorithm is as follows. 



Algorithm Alternating Split AL 
1. Set k = 0, choose fi > 0, uo, vo, and do. 
repeat 

Ufc+i £ argmin u /i(u) + f ||u - v fc 



dfc||l 



Vfc + i £ argmin v / 2 (v) + f ||u fc+ i - v - d» 
dfc+i <— dfc — Ufc + i + Vfc + i 
k <- k + 1 
until stopping criterion is satisfied. 



Problem (1) has the form (3) where /i is quadratic, thus Step 3 
consist in solving a linear system of equations. We will return to the 
particular form of this system in the next section. With / 2 = r<f>, a 
regularizer, Step 4 corresponds to applying a shrinkage/thresholding 
function, that is, Vfc + i = 4 , T 0/ M (uk+i — dfc) , usually a computa- 
tionally inexpensive operation. 



3. PROPOSED METHOD 

3.1. Constrained Optimization Formulation and Algorithm 

Performing the VS on problem (1) yields the following constrained 
formulation: 

mm !||HW/3-y||! + r0(0) 

P,8 (10) 

s.t. (3 = 9. 

This VS decouples the quadratic non-separable term ||HW/3 — y 
from the non-quadratic term 4>(6), to deal with the non-separability 
of the quadratic data term. In contrast, split-Bregman methods use a 
splitting to avoid non-separability of the regularizes 

Problem (10) has the form (4), with u = /3, v = 0, /i(u) = 
(l/2)||HW/3-y|||,and/ 2 (v) =rcj>{0). Applying this translation 
table to the Alternating Split AL algorithm presented Section 2.3, we 
obtain the following algorithm. 

Algorithm Split AL Shrinkage Algorithm 

1. Set k = 0, choose fi > 0, /3 , Oo, and do. 

2. repeat 

3. p' k = e k + d k 

4. k+1 eargrmn||HW/3-y||l+ M ||/3-/3' fc ||l 

5. 0' k = p k+1 - d k ^ 

6. k +i = 

7. d k+1 <— d k — (3 k+1 + Ok+i 

8. k^k + 1 

9. until stopping criterion is satisfied. 

Since Step 4 is a strictly convex quadratic problem, its solution 
is unique and given by 

0k+i = (W T H T HW + M l) _1 (w T -a T y + p(3' k ) . (11) 

In the next subsection, we show how f3 k+1 can be efficiently com- 
puted. Note that (W T H T HW+/il) is a regularized (by the addi- 
tion of pi) version of the Hessian of | ||HW/3 — y 111- 

3.2. Computing f3 k+1 

Assume that W is a normalized tight (Parseval) frame, i.e. , W W T = 
I (although possibly W T W 7^ I), and that H is the matrix represen- 
tation of a convolution, i.e., products by H or H T can be computed 
in the Fourier domain, with 0(n log n) cost via the FFT 

The assumptions in the previous paragraph will enable us to 
compute the matrix inversion in (1 1), even if it is not feasible to ex- 
plicitly form matrix HW. Using the Sherman-Morrison- Woodbury 
inversion formula, (11) becomes 

(3 k+1 = ^I-WVlHWWV+plj^Hwjrt 

= ^I-W T H T (HH T + filj _1 Hwjr t (12) 

where Vk = (W T H T y + p (3' k ). Furthermore, since H is the ma- 
trix representation of a convolution, (12) can be written as 

/3 fc+ i = J (l-W T Fw)r fc , (13) 

where 

F = XJ H H* (|D| 2 DU, (14) 



U and U are the matrix representations of the forward and in- 
verse discrete Fourier transform (DFT), and D is a diagonal matrix 
containing the DFT of the convolution represented by H. Notice that 
the product by F corresponds to applying a filter in the DFT domain, 
which can be done using FFT algorithms with 0(n log n) cost. No- 
tice also that H T W T y can be precomputed. When the products by 
W T and W are direct and inverse tight frame transforms for which 
fast algorithms exist, the leading cost of each application of (13) will 
be either 0(n log n) or the cost of these frame transforms (usually 
also 0(n log n)). 

Finally, the complete algorithm, which we term SALSA (split 
augmented Lagrangian shrinkage algorithm) is as follows. 

Algorithm SALSA 

1. Initialization: set k — 0; choose p > 0, (3 , do, do; 

2. compute y = W T H T y 

3. compute F = U ff D* (|D| 2 + ^1) _1 DU 

4. repeat 

5. p' k <-e k + d k 

6. r k <- y + pj3' k 

7- P k+1 «- £ (I - W T F fc W) r fc 

8. 6' k «- f3 k+1 - d k 

9. e k+1 «- * rHll {# k ) 

10. d fe+ i <— d fc - (3 k+1 + 6 k+1 

11. k^k + 1 

12. until stopping criterion is satisfied. 



4. EXPERIMENTS 

We consider five standard image deconvolution benchmark problems 
[13], summarized in Table 1, all on the well-known Cameraman im- 
age. The regularizer is (f>((3) = ||/3||i, thus is a soft threshold. 
In all the experiments, W is a redundant Haar wavelet frame, with 4 
levels, and the blur operator H is applied via the FFT. The regular- 
ization parameter r in each case was hand tuned for best improve- 
ment in SNR. The value of p, for fastest convergence was found to 
differ in each case, but a good rule of thumb, used in all the experi- 
ments, is p = O.lr. We compare SALSA with current state of the 
art methods: TwIST [3], SpaRSA [22], and FISTA [2], in terms of 
the time taken to reach the same value of the objective function. Ta- 
ble 2 shows the CPU times taken by each of the algorithms in each of 
the experiments. Figure 4 shows the plots of the objective function 
i 1 1 HW/3 — y 1 1 + t 1 1 P 1 1 1 , evolving over time, in experiments 1 , 2B , 
and3A. 



Table 1. Details of the image deconvolution experiments. 



Experiment 


blur kernel 


a 2 


1 


9x9 uniform 


0.56^ 


2A 


Gaussian 


2 


2B 


Gaussian 


8 


3A 
3B 


hi S = i/(i + ; 2 +j 2 ) 
hij = 1/(1 + i 2 +f) 


2 
8 



5. CONCLUSIONS 

We have proposed a fast algorithm for frame-based image deconvo- 
lution, based on variable splitting and solving the constrained opti- 
mization problem through an augmented Lagrangian scheme. Ex- 
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(a) (b) (c) 

Fig. 1. Objective function evolution: (a) 9 x 9 uniform blur, a = 0.56; (b) Gaussian blur, a 2 = 8; (c) hij — 1/(1 + i 2 + j 2 ) blur, a 2 = 2. 



Table 2. CPU times (in seconds) for the various algorithms. 



Experiment 


TwIST 


SpARSA 


FISTA 


SALSA 


1 


50.2969 


42.0469 


64.2344 


4.000 


2A 


30.7656 


40.6094 


61.7031 


4.03125 


2B 


14.4063 


6.92188 


15.0781 


1.9375 


3A 


23.5313 


17.0156 


33.7969 


2.60938 


3B 


8.1875 


6.17188 


18.0781 


1.89063 



perimental results with l\ regularization show that our new algo- 
rithm outperforms existing state-of-the-art methods in terms of com- 
putation time, by a considerable factor. Future work includes the ap- 
plication of SALSA to other inverse problems, namely compressed 
sensing and reconstruction with missing samples. 
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